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The displacement of star images by atmospheric refraction observed by an Earth-bound telescope 
is dominated by a familiar term proportional to the tangent of the zenith angle and proportional to 
the refractivity at the ground. 

The manuscript focuses on the torsion of the ray path through the atmosphere in a model of 
ellipsoidal atmospheric layers above the Earth surface, induced by the two slightly different principal 
curvatures along N-S and E-W pointing directions, depending on the geodetic latitude of the 
telescope site. This symmetry breaking effects apparent places in the sub-milliarcsecond range 
at optical and infrared wavelengths. 

PACS numbers: 95.75.Mn, 95.10.Jk 

Keywords: Atmospheric Refraction, Astronomy, Apparent Place, Prolate Ellipsoid 



Or 
i 

O 



> 
o 



o 

OS 

o 



X 



I. RAY PATHS FROM SNELL'S LAW 

Tracing of rays of stellar light through the Earth atmo- 
sphere is a repeated application of Snell's law of refrac- 
tion, as the refractive index n increases from n — 1 = 
high above the atmosphere to n — 1 of the order of 
2 x 10~ 4 — perhaps 3 x 10~ 4 at sea level — at the tele- 
scope site. The angle between their direction t at ar- 
rival at the telescope and the topocentric zenith is smaller 
than in a calculation with vacuum replacing the atmo- 
sphere (E Bi IS EH E3] ; the atmosphere is a gradient lens 
BUS El- The ray path is curved. Snell's law of the 
relation between refractive index n and incidence angle 
if) relative to the normal to the layered atmosphere is 



nsin?/> = const. 
The first differential of this equation is 

An sin ip + n cos tpAtj) = 0, 
or upon division through cos ip 

An tan tp = —nArp, 



(1) 



(2) 



(3) 



which explains in conjunction with n — 1 <C 1 why the 
difference between true and apparent positions, the ac- 
cumulation of all Aip, basically turns out proportional to 
the tangent of the zenith angle. 

The continuum version is equivalent to solving the dif- 
ferential equation which couples the gradient of the re- 
fractive index to the ray's curvature @, chap. 3.2] 



d ( dr\ ^ 

as \ as . 



(4) 



for the path r as a function of path length s. This can be 
transformed to general curvilinear coordinates and 



has been made explicit for the Earth ellipsoid by Yat- 
senko The present work demonstrates that the as- 
phericity changes apparent places by less than one milli- 
arcsecond, such that a simple finite-element computation 
suffices for all but some ambitious astrometry projects. 



II. NUMERICAL SIMULATION 

A. Oblate Ellipsoid Coordinates 

The Cartesian coordinates of a position r are related to 
the geodetic latitude c/>, geodetic longitude A and height 
h above the ellipsoid as [8| 



r = 




where 



N(cj>) 



h) cos 4> cos A 
h) cos (j> sin A 
- e 2 ) + h\ sin0 



yl — e 2 sin 2 < 



(5) 



(6) 



is the distance to the Earth axis along the surface normal. 
p e ~ 6378 km and e sa 0.0818 arc equatorial radius and 
eccentricity |13| . 



B. Refractive Index of Altitude 

The simulation has been performed with an exponen- 
tial model of the air susceptibility x(h) as a function of 
altitude h, which changes in discrete steps at altitudes 
hi such that the kinks n(/ij+i) — n(hi) are approximately 
equal across each interface: 



n (h) = y/l + Xu hi < h < hi 



(J) 



*URL: http://www.strw.leidenuniv.nl/~mathar; Electronic ad- 
dress: mathar@strw.leidenuniv.nl 



This is achieved by defining a set of heights loga- 
rithmic function of index i, using each second to specify 
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interface locations hi, and using the complementary set 
of interlaced heights to specify the x%- 

a, = K\og(-^—), Q<i<ioo] (8) 

K = d2i-i] 1 < i < W 2 ; (9) 
X* = Xoe- a2 ' /K , 0<z< W2;XW2 = °- ( 10 ) 

The scale height was set to K = 9.6 km and the ground 
layer value to xo — 4 x 10~ 4 . The results are practicably 
stable once the number of layers, that is ioo/2, is set to 
values of order 10 or 20. 



C. Finite Element Segmentation 

Rays are traced backwards from the position ro of the 
telescope to the vacuum above the troposphere. This 
removes the need to implement some adaptive forward 
searching techniques for the alternative handling which 
would trace the ray in its natural direction. 

The implementation is much simpler than writing 
down the associated first order differential equation of 
the eikonal in the curvilinear coordinates @. 

At each layer interface, a topocentric coordinate sys- 
tem is spanned by the unit vectors e^, and [l2| 

d d d 

The positions along the straight paths inside the layers 
are defined as 

r = n + tn (12) 

via a direction r and distance parameter t. Given the al- 
titude hi + \ of the next layer as defined by the layer model 
above, the unknowns t, <f>i+i and A»+i which constitute 
rj + i are calculated by solving 

/ [N(<f> i+ i) + hi+i] cos0 l+ i cos Ai+i \ 
ri+<x i= [N((j) l+1 ) + h i+ i] cos^+i sinA i+ i .(13) 
\ [JV(0 i+ i)(l - e 2 ) + h i+1 ] Bin^+i / 

This nonlinear system of 3 equations and 3 unknowns is 
reduced to 2 equations and 2 unknowns by taking the 
sum of the squares of the first two components, which 
eliminates Aj+i, and then solved by a Newton method. 

Once the ftf+i, 4>i+i an d A^+i are known, a local triad 
(ej,,e\,eh) is set up at the interface point Vi+\. The 
direction Tj is virtually decomposed as 

Tj = ae^i+i + (3e x ,i + i + je h ,i+i (14) 

that is, the angle of incidence ^ is computed as the pro- 
jection 

n ■ e/i, i+ i = cos^lTiHe^^i+il = 7|e M+ i| 2 . (15) 



Snell's law fixes the exit angle x 

n i+ i sinx = m sinip, (16) 
and the new direction tj+i is constructed by switching 
in a new projection 7' for the old 7 as 

7 tan?/; = 7'tanx; (17) 
t 1+ i = ae^ i+ i+Pex,i+i+ r /eh,i+i. (18) 

The procedure is iterated with the parametrization r,; + i + 
tTi + \ for the path inside the next layer. The decompo- 
sition (|14p is virtual in the sense that a, e^j+i, j3 and 
©A.j+i are not actually computed; the new direction is an 
update of the old by adding the Cartesian components 
of (7' — "/)eh,i+i to the Cartesian components of 73. 

The format (fl~4"l) defines the apparent zenith angle zq 
and apparent azimuth Aq in the ground layer 

to ■ e^ = sinz o cosv4 o |e0 ; o| 2 ; (19) 
To-e A ,o = -sinz sm,4o|eA,o| 2 ; (20) 
to • e/,,0 = coszole/^ol 2 , (21) 

where the azimuth A is counted North-over- West. The 
simulation produces the true zenith angle z and true az- 
imuth A by the equivalent projection of t^, the direc- 
tion at exit above the troposphere, onto the axes of the 



topocentric system at the telescope, i.e., 

Too • e^o = sinzcos^|e0 :O | 2 ; (22) 

Too-eA,o = -sinzsinA|e A ,o| 2 ; (23) 

Too ■ e/,,0 = cosz|e M | 2 . (24) 



III. RESULTS 
A. Apparent Azimuth 

The difference z — zq is the familiar change in zenith 
angle. The non-zero difference A — Aq is a measure of the 
ray torsion Q , which vanishes in the limit of zero eccen- 
tricity. This change of apparent azimuth is illustrated 
in Fig. Q] for the Earth parameters quoted above. The 
zenith angle zo and azimuth Aq are coded in circular co- 
ordinates at the base of each sub-plot: in the center, zq is 
zero, and it grows outwards in steps of 5° such that the 
maximum zq = 60° is reached at the rim of the graph. 
The azimuth directions N and W are indicated by arrows. 

The presence of torsion such that the rays do not stay 
in a single plane is a side effect of the two different cur- 
vatures along the N-S and E-W directions set up by this 
atmospheric model (App.[A|. The curvature is larger for 
pointing into N-S than for pointing into E-W directions. 
Larger curvature implies smaller air mass, i.e., smaller 
angles of incidence relative to the normals of the air lay- 
ers, and implies smaller optical path length. 
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FIG. 1: The difference A — Ao between true and apparent azimuth angles in units of milli-arcseconds for eight different geodetic 
latitudes <f> of the telescope as a function of pointing direction (zenith angle zero to 60° growing outwards; azimuth with a 
deliberate gap in the South direction at ±180°). 
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(f>=10° <t>=20° 




FIG. 2: The difference in star zenith angles obtained with a prolate ellipsoid compared with the computation with a spherical 
Earth, (z — zo)\ e — (z — zo)\ e =o, m units of milli-arcseconds, for eight geodetic latitudes <j) plotted as a function of zo < 60°, 
-175° < Ao < 175°. 
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Since Snell's law is a result of the eikonal minimization 
(Fermat's principle) which lets light select minimum op- 
tical path length, the net effect is to steer light closer to 
the N if it has a northern component, and closer to S if it 
has southern component. In our convention — azimuths 
wrapped around such that \A\ < 180° — , the directions 
of set azimuths \ A\ < A n if \ A \ < 90°, and azimuths 
\A\ > A a if |A | > 90°. This explains why A - A a < 
for NW or SE directions, and why A — A > for SW 
or NE directions, as shown in Fig. [T] The effect vanishes 
approaching the geodetic poles, <j> — > ±90°, because the 
two curvatures equalize there (Eq. (jAlOj) ). 

The discontinuity of the plots in Fig. Q] in the centers 
near z = is the standard artifact of polar coordinates. 
The same reasoning applies which leads to the multipli- 
cation of right ascension components of proper motions 
by the cosine of the declinations. 



IV. SUMMARY 

Ray paths through the atmosphere in the geometric op- 
tics approximation change their direction in zenith angle 
and azimuth if the refractive indices in the atmosphere 
layers are shells modeled with a symmetry of the prolate 
ellipsoid of the Earth. The change in azimuth direction 
is less than one milli-arcsecond for pointing less than 60° 
away from the zenith at typical optical wavelengths, com- 
parable to the diffraction limit of a 170 m telescope, and 
therefore much smaller than the standard refractivity ef- 
fect (and others not discussed here [l5j]) which changes 
apparent positions on the scale of 60 arcseconds under 
the same conditions. 



APPENDIX A: CURVATURES 



B. Apparent Zenith 

In Fig. [2] we compare the change in zenith angles of 
Too as the ellipticity e is switched off to e = 0, keep- 
ing p e the same. This effect is typically buried under 
the standard effect up to 60 arcseconds. We see that 
for telescope sites near the equator, small |</>|, the zenith 
angles for pointing in E or W direction are the same as 
calculated for the spherical model, whereas the N or S di- 
rections have smaller refractivities. Fig.[2]just tracks the 
changes in the two principal curvatures (|A10[) and (|A12[) 
and finishes with positive values near \<fi\ — > 90° because 
the curvature of the prolate model near the poles, radius 
p p > p e , is smaller than the curvature of the spherical 
model kept at p e . 

For (j> < 54° [up to where the line of K2(4>) meets the 
horizontal defined by Ki(0) in Fig. [3], one can find az- 
imuths A for which the difference vanishes as the cur- 
vature k = sin 2 Ak% + cos 2 Ak,2 of Euler's formula equals 
the equatorial curvature. This explains the saddle shapes 
in Fig. [2] while <j> < 54°. Evidently, the figure would look 
different if the reference were an Earth sphere with a ra- 
dius selected from a mean (Gauss) curvature depending 
on (p, for example. 

This effect results from any accurate calculation of the 
standard zenithal apparent places with variable earth ra- 
dius. The amount of 1 milli-arcsecond for <j> = 80° in 
Fig. [2] is expected from my earlier calculation [ill, Fig. 
4], which stated that the effect of replacing the zero 
curvature model by a model with the standard curva- 
ture is 0.5% of z — zq, which accumulates 0.3 arcseconds 
supposed z — zq ~ 60 arcseconds. Because the factor 
l/\/l — e 2 of (jAlOjl introduces a relative difference of 
0.3% between equatorial and polar radii or curvatures, 
linear extrapolation predicts indeed a change of up to 
another 1 milli-arcsecond in apparent zenith angle if tele- 
scopes are moved between equator and poles [l8j. 



The Fundamental Parameters of the first form of the 
surface of constant altitude h for prolate ellipsoids are 

EI 



E = {N + hf cos 2 , 
F = 0; 

G = (M + hf. 



(Al) 
(A2) 
(A3) 



The Fundamental Parameters of the second form of the 
surface are 



L = -{N + h)cos 2 <f>; 
M = 0; 

N = -(M + h), 



where 



M(4>) = N(4>) 



1 — e 2 sin 2 i 



(A4) 
(A5) 
(A6) 



(A7) 



The six Fundamental Parameters are adorned with an 
overbar to set and M apart from the distances N and 
M in ([6]) and (|A7[) . The two principal curvatures Ki,2(<A) 
are the roots of the quadratic equation 



(EG - F z )k z - (EN - 2FM + GL)n + (LN - M z ) = 0. 

(A8) 

Insertion of the 6 parameters yields 

(N + h)(M + K)k 2 + (M + N + 2K)k +1 = 0. (A9) 

The two solutions K\ = —l/(N+h) and K2 = — 1/(M+ 
h) are shown in Fig. [3J In particular: 



• Above the poles, s'm<fi — ±1, N — M and 

1 



fSi, 2 (0 = ±7r/2) = 



(A10) 



so the curvatures are derived from the scmiminor 
axis p p — p e /Vl — e 2 [i"4j . 
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FIG. 3: The principal curvatures from (|A9|I at h = as a 
function of geodetic latitude. 



Above the equator, sin<^> = 0, N = p e , M — p e (l 
e 2 ) and 



«ifa = 0) = - 
«a(^ = 0) : 



Pe + h 1 



p e (l-e 2 )+h' 



(All) 
(A12) 



which associates with the equatorial radius of cur- 
vature p e and the polar radius of curvature p e (l — 
e 2 ) 0. 



In the main text we call \k\ the curvature; coherent 
with the standard nomenclature for closed shell surfaces, 
the negative sign is maintained in this appendix. 
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